###Space, Time, Space vs. Time#####
###02/01/2022######################
###Generate Data##################

#Generic working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#Get packages, attach them
pks=c("spdep", "lmtest","Matrix","spatialreg","tidyverse")
needed=setdiff(pks,installed.packages())

if( length(needed) > 0){
  options(Ncpus=parallel::detectCores() %/% 1.25)
  install.packages(needed)
}

sapply(pks, library, character.only=T)


options(scipen=999)
set.seed(5862685)

N = 50
t = 20  

tr = 250


#Generate W matrix
x.coord = runif(N,0,100)
y.coord = runif(N,0,100)

points = cbind(x.coord,y.coord)

dnb = knearneigh(points, 5, longlat = NULL)
dsts = knn2nb(dnb)
lw = nb2listw(dsts, style="W", zero.policy=TRUE)
W = listw2mat(lw)


#Pre-alocate for storage
dp = 7 
beta.results.static = matrix(NA, nrow=5, ncol=dp)
beta.results.ldv = matrix(NA, nrow=5, ncol=dp)
beta.results.sar = matrix(NA, nrow=5, ncol=dp)
beta.results.stadl = matrix(NA, nrow=5, ncol=dp)
phi.results.ldv = matrix(NA, nrow=5, ncol=dp)
phi.results.stadl = matrix(NA, nrow=5, ncol=dp)
rho.results.sar = matrix(NA, nrow=5, ncol=dp)
rho.results.stadl = matrix(NA, nrow=5, ncol=dp)

lm.results.ldv = matrix(NA, nrow=5, ncol=dp)
lm.test = matrix(NA, nrow=1, ncol=dp)

#Fixed elements 
beta = 2
rho_x = 0.3
phi_x = 0.6 

X = rnorm(N)
I = diag(N)
inv_x = solve(I - (rho_x * W)) 
X_0 = inv_x%*%(X)

IT = diag(N*t)
IT_1 = diag(N*(t-1)) 
TL = cbind(rbind((matrix(0, N, N*(t-1))),IT_1),(matrix(0, N*t, N)))

IT_1_resid = diag(N*(t-2)) 
TL_resid = cbind(rbind((matrix(0, N, N*(t-2))),IT_1_resid),(matrix(0, N*(t-1), N)))

X_0 = rep(X_0,t) 
WT = as.matrix(bdiag(replicate(t,W,simplify=F)))

x_mult = solve(IT - phi_x*TL - rho_x*WT) 
X = x_mult%*%(X_0 + rnorm(N*t)) 

XTL = TL%*%X
WX = WT%*%X

WT_1 = as.matrix(bdiag(replicate(t-1,W,simplify=F)))
wtl = mat2listw(WT_1, row.names = NULL, style="W")

#Stochastic elements
for(k in 1:6){
for(r in 1:dp){  

rho_y = (r-1)*0.05
theta = 0

phi_y = (k-1)*0.10
gamma = 0
tim_spa = solve(IT - phi_y*TL - rho_y*WT)



beta.coef = matrix(NA, nrow=tr, ncol=4)
beta.se = matrix(NA, nrow=tr, ncol=4)
phi.coef = matrix(NA, nrow=tr, ncol=2)
phi.se = matrix(NA, nrow=tr, ncol=2)
rho.coef = matrix(NA, nrow=tr, ncol=2)
rho.se = matrix(NA, nrow=tr, ncol=2)
res.coef = matrix(NA, nrow=tr, ncol=1)
res.se = matrix(NA, nrow=tr, ncol=1)
test = matrix(NA, nrow=tr, ncol=1)

for(i in 1:tr){  
  Y = tim_spa%*%(beta*X +  rnorm(N*t))
  
  LDV = TL%*%Y
  LDV[LDV == 0] <- NA
  
  
  id = rep(seq(1,N,1),t) 
  time = rep(1:t, each=N)
  
  DGP = as.data.frame(cbind(Y,X,LDV,WX,XTL,id,time))
  DGP = na.omit(DGP)
  
  colnames(DGP) <- c("Y","X","LDV","WX","XTL","ID","TIME")
  

  
  stat = lm(Y~X) #STATIC 
  ldv = lm(Y~X+LDV, data=DGP) #LDV
  sar = lagsarlm(Y~X,data=DGP, wtl) #SAR 
  stadl = lagsarlm(Y~X+LDV,data=DGP, wtl) #STADL(1,1)
  
  res = residuals(ldv)
  res1 = TL_resid%*%res
  LM = lm(res~X+LDV+res1, data=DGP)
  test[i,] = coef(summary(LM))[4,4]
  

  beta.coef[i,] = c(coef(summary(stat))[2,1],coef(summary(ldv))[2,1],summary(sar)$coefficients[2], summary(stadl)$coefficients[2]) 
  beta.se[i,] = c(coef(summary(stat))[2,2],coef(summary(ldv))[2,2],summary(sar)$rest.se[2], summary(stadl)$rest.se[2]) 
  
  
  phi.coef[i,] = c(coef(summary(ldv))[3,1], summary(stadl)$coefficients[3])
  phi.se[i,] = c(coef(summary(ldv))[3,2], summary(stadl)$rest.se[3])
  
  rho.coef[i,] = c(summary(sar)$rho, summary(stadl)$rho)
  rho.se[i,] = c(summary(sar)$rho.se, summary(stadl)$rho.se)

  res.coef[i,] = coef(summary(LM))[4,1]
  res.se[i,] = coef(summary(LM))[4,2]
}

beta.bias = colMeans(beta.coef, na.rm = FALSE, dims = 1) - beta 
beta.ase = colMeans(beta.se, na.rm = FALSE, dims = 1)
beta.std = apply(beta.coef, 2, sd)
beta.rmse = matrix(NA, nrow=1, ncol=4)
beta.cp = matrix(NA, nrow=1, ncol=4)
for(m in 1:4){ 
  beta.rmse[m] = (sum((beta.coef[,m] - beta)^2))/tr
  beta.UB = beta.coef[,m]+1.96*beta.se[,m]
  beta.LB = beta.coef[,m]-1.96*beta.se[,m]
  beta.cp[m] = sum(ifelse((beta>= beta.LB)& (beta<= beta.UB),1,0)) 
}

beta.results.static[,r] = rbind(beta.bias[1],beta.ase[1],beta.std[1],beta.rmse[1],beta.cp[1])
beta.results.ldv[,r] = rbind(beta.bias[2],beta.ase[2],beta.std[2],beta.rmse[2],beta.cp[2])
beta.results.sar[,r] = rbind(beta.bias[3],beta.ase[3],beta.std[3],beta.rmse[3],beta.cp[3])
beta.results.stadl[,r] = rbind(beta.bias[4],beta.ase[4],beta.std[4],beta.rmse[4],beta.cp[4])


phi.bias = colMeans(phi.coef, na.rm = FALSE, dims = 1) - phi_y 
phi.ase = colMeans(phi.se, na.rm = FALSE, dims = 1)
phi.std = apply(phi.coef, 2, sd)
phi.rmse = matrix(NA, nrow=1, ncol=2)
phi.cp = matrix(NA, nrow=1, ncol=2)
for(m in 1:2){ 
  phi.rmse[m] = (sum((phi.coef[,m] - phi_y)^2))/tr
  phi.UB = phi.coef[,m]+1.96*phi.se[,m]
  phi.LB = phi.coef[,m]-1.96*phi.se[,m]
  phi.cp[m] = sum(ifelse((phi_y>= phi.LB)& (phi_y<= phi.UB),1,0)) 
}
phi.results.ldv[,r] = rbind(phi.bias[1],phi.ase[1],phi.std[1],phi.rmse[1],phi.cp[1])
phi.results.stadl[,r] = rbind(phi.bias[2],phi.ase[2],phi.std[2],phi.rmse[2],phi.cp[2])


lm.test[,r] = sum(ifelse(test<0.05,1,0))

rho.bias = colMeans(rho.coef, na.rm = FALSE, dims = 1) - rho_y
rho.ase = colMeans(rho.se, na.rm = FALSE, dims = 1)
rho.std = apply(rho.coef, 2, sd) 
rho.rmse = matrix(NA, nrow=1, ncol=2)
rho.cp = matrix(NA, nrow=1, ncol=2)
for(m in 1:2){ 
rho.rmse[m] = (sum((rho.coef[,m]- rho_y)^2))/tr
rho.UB = rho.coef[,m]+1.96*rho.se[,m]
rho.LB = rho.coef[,m]-1.96*rho.se[,m]
rho.cp[m] = sum(ifelse((rho_y>= rho.LB)& (rho_y<= rho.UB),1,0)) 
}
rho.results.sar[,r] = rbind(rho.bias[1],rho.ase[1],rho.std[1],rho.rmse[1],rho.cp[1])
rho.results.stadl[,r] = rbind(rho.bias[2],rho.ase[2],rho.std[2],rho.rmse[2],rho.cp[2])

}

filename <- paste("FINALphi", phi_y, sep = "_")
filename <- paste(filename, ".Rdata", sep = "")
save.image(file=filename)
}

